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Real-life  quantum  computers  are  inevitably  affected  by  intrinsic  noise  resulting  in  dissipative  nonunitary 
dynamics  realized  by  these  devices.  We  consider  an  open-system  quantum  annealing  algorithm  optimized 
for  such  a  realistic  analog  quantum  device  which  takes  advantage  of  noise-induced  thermalization  and 
relies  on  incoherent  quantum  tunneling  at  finite  temperature.  We  theoretically  analyze  the  performance  of 
this  algorithm  considering  a  p-spin  model  that  allows  for  a  mean-freld  quasiclassical  solution  and,  at  the 
same  time,  demonstrates  the  first-order  phase  transition  and  exponential  degeneracy  of  states,  typical 
characteristics  of  spin  glasses.  We  demonstrate  that  finite-temperature  effects  introduced  by  the  noise  are 
particularly  important  for  the  dynamics  in  the  presence  of  the  exponential  degeneracy  of  metastable  states. 

We  determine  the  optimal  regime  of  the  open-system  quantum  annealing  algorithm  for  this  model  and  find 
that  it  can  outperform  simulated  annealing  in  a  range  of  parameters.  Large-scale  multiqubit  quantum 
tunneling  is  instrumental  for  the  quantum  speedup  in  this  model,  which  is  possible  because  of  the  unusual 
nonmonotonous  temperature  dependence  of  the  quantum-tunneling  action  in  this  model,  where  the  most 
efficient  transition  rate  corresponds  to  zero  temperature.  This  model  calculation  is  the  first  analytically 
tractable  example  where  open-system  quantum  annealing  algorithm  outperforms  simulated  annealing, 
which  can,  in  principle,  be  realized  using  an  analog  quantum  computer. 
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I.  INTRODUCTION 

Quantum  computing  hardware  is  affected  by  a  substantial 
level  of  intrinsic  noise  and  therefore  naturally  realizes  dis¬ 
sipative  quantum  dynamics  [1,2].  Optimization  algorithms, 
where  a  configuration  of  a  binary  string  x  minimizing  a  given 
(energy  or  cost)  function  f(x)  is  sought  for,  naturally  extract  a 
computational  advantage  from  the  irreversible  dissipative 
dynamics  and  could  therefore  be  readily  implemented  on  a 
number  of  existing  hardware  platforms  [3,4].  More  specifi¬ 
cally,  quantum  annealing  [5]  (QA)  is  a  quantum  analog  of  the 
widely  applied  classical  simulated  annealing  algorithm  [6] 
(SA),  a  heuristic  solver  of  NP-hard  (non-deterministic  poly¬ 
nomial-time  hard  [7])  optimization  problems  [5,8-11],  with 
quantum  fluctuations  playing  the  role  analogous  to  thermal 
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fluctuations  in  simulated  annealing.  NP-hard  optimization 
problems,  such  as  finding  a  ground-state  spin  configuration  of 
a  spin  glass,  are  often  characterized  by  an  energy  landscape 
with  a  large  number  of  local  minima  separated  by  extensive 
energy  barriers.  Dissipative  dynamics  realized  by  the  open- 
system  quantum  annealing  provides  an  efficient  mechanism 
for  thermalization  within  domains  of  attraction  of  local 
minima.  For  an  efficient  search  of  the  configuration  space, 
the  barriers  separating  different  domains  of  attraction  have  to 
be  overcome,  which  may  proceed  via  thermal  excitation  or  a 
quantum-tunneling  process.  The  performance  of  the  open- 
system  quantum  annealing  algorithm  is  therefore  character¬ 
ized  by  a  set  of  relaxation  rates  associated  with  such  processes, 
as  opposed  to,  for  example,  the  spectral  gaps,  as  is  the  case  for 
an  adiabatic  quantum  algorithm  [10,12,13]. 

The  longest  relaxation  times  correspond  to  the  often 
exponentially  slow  transitions  between  local  minima  sep¬ 
arated  by  extensive  potential  banners.  Unitary  dynamics  of 
a  pair  of  such  states  conesponds  to  the  switching  rate  of  the 
order  of  the  matrix  element  A,  which  in  the  presence  of  an 
extensive  barrier  may  scale  exponentially  with  the  system 
size  A  oc  exp(— const  x  N )  (here,  N  is  the  number  of  qubits 
in  the  system).  However,  fast  dissipative  relaxation  within  a 
domain  of  attraction  of  a  local  minima  due  to  the  hardware 
noise  introduces  a  lifetime  or  level  width  W.  This  fast  local 
relaxation  strongly  suppresses  the  coherent  superposition 
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of  the  states  localized  in  different  local  minima  when 
W  :»  A.  Nevertheless,  the  incoherent  quantum  tunneling 
is  possible  in  the  presence  of  such  strong  dissipation,  where 
the  transition  rate  is  described  by  the  Fermi  golden-rule-type 
expression  a  A2  «:  A.  This  is  the  regime  likely  realized  in  a 
large-scale  quantum  annealer  [3,14].  It  is  an  open  question 
whether  such  incoherent  extensive  quantum  tunneling  may 
provide  a  more  efficient  mechanism  for  searching  the 
configuration  space  as  compared  to  classical  simulated 
annealing  relying  on  thermal  excitation.  Initial  numerical 
analysis  of  the  two-dimensional  (2D)  spin-glass  problem 
[15]  suggested  a  superior  scaling  of  the  outcome  of  finite- 
temperature  quantum  annealing  compared  to  simulated 
annealing  [11].  However,  this  observation  turned  out  to  be 
an  artifact  of  the  numerical  discretization  scheme  and 
therefore  cannot  be  reproduced  using  analog  hardware  [17]. 

In  this  paper,  we  provide  an  analytically  tractable  example 
where  the  speedup  originates  from  incoherent  quantum 
tunneling.  We  consider  a  system  of  Ising  spins  interacting 
each  with  each  other  with  a  p-body  interaction  of  equal 
strength,  a  model  often  referred  to  as  a  p-spin  model.  This 
model  allows  a  quasiclassical  Wentzel-Kramers-Brillouin 
(WKB)  description  [18,19],  where  the  expansion  is  per¬ 
formed  in  1  /N  rather  than  the  more  usual  h,  and,  at  the  same 
time,  it  demonstrates  key  features  characteristic  of  a  range  of 
complex  (NP-hard)  optimization  problems,  such  as  the  first- 
order  phase  transition  (for  p  >  3)  and  an  exponentially  small 
gap  between  the  ground  and  excited  states  [20].  Crucially, 
the  metastable  state  realized  in  this  model  is  characterized 
by  an  exponential  degeneracy,  whereas  the  ground  state  is 
unique.  Such  entropic  imbalance  is,  in  fact,  typical  for 
low-energy  states  in  the  spin-glass  phase,  and  it  strongly 
affects  the  low-temperature  dynamics  of  the  system  in  both 
quantum  and  classical  cases.  The  effect  of  entropic  imbal¬ 
ance  is  the  main  focus  of  our  analysis  in  this  paper. 

We  show  that  the  scaling  of  the  optimal  QA  computation 
time  (allowing  for  repeated  runs  of  the  algorithm)  is 
determined  by  the  quantum-tunneling  amplitude  at  a  single 
point  in  the  algorithm,  the  so-called  freezing  point,  after 
which  quantum  (or  thermal)  fluctuations  are  relatively  weak 
and  the  transitions  over  or  through  the  barrier  are  no  longer 
likely.  We  find  that,  because  of  the  effect  of  the  entropy 
associated  with  the  metastable  state,  the  optimal  quantum¬ 
tunneling  rate  is  achieved  at  vanishing  temperature;  i.e., 
raising  the  temperature  may  reduce  the  quantum-tunneling 
rate.  This  is  in  contrast  to  the  usual  intuition  about  a 
quantum-mechanical  particle  trapped  in  a  nondegenerate 
metastable  potential  well,  where  the  escape  rate  monoto¬ 
nously  increases  with  temperature.  The  optimal  QA  regime 
corresponds  to  the  fastest  tunneling  rate  and  therefore 
vanishing  temperature  as  well.  Comparing  the  optimal 
computation  time  of  QA  obtained  in  this  regime  with  that 
of  S  A  for  a  range  of  the  potential  banner  shapes,  we  find  that 
QA  could  outperform  SA  under  certain  circumstances,  thus 
providing  a  polynomial  (rather  than  exponential)  speedup. 


The  physical  mechanism  of  the  quantum  speedup  that 
we  find  is  distinct  from  the  usual  intuition  of  quantum 
fluctuations  overcoming  thin  and  tall  banners  more  effi¬ 
ciently.  It  is  a  generic  mechanism,  and  we  expect  it  to 
manifest  in  more  complex  quantum  models  where  quantum 
fluctuations  are  introduced  by  transverse  spin-spin  or 
multispin  interaction  terms  [21-23]  in  addition  to  the 
transverse  field.  We  emphasize  that  the  speedup  mechanism 
we  find  is  not  limited  to  the  models  that  can  be  efficiently 
simulated  with  quantum  Monte  Carlo  dynamics  on 
classical  computers  [24,25]  (such  as  the  transverse-field 
Ising  model).  The  regime  considered  here  could,  in 
principle,  be  reproduced  using  an  analog  quantum  annealer. 

Before  outlining  the  formal  calculation,  we  discuss  the 
qualitative  picture  of  the  effect  of  the  entropy  imbalance 
between  the  ground  and  metastable  states  on  the  efficiency  of 
simulated  and  quantum  annealing.  In  the  presence  of  the 
entropy  imbalance,  SA  computation  time  scales  exponen¬ 
tially  with  the  system  size  N.  This  can  be  understood 
intuitively  by  considering  the  performance  of  SA  applied 
to  a  model  demonstrating  a  first-order  phase  transition  into 
a  state  characterized  by  an  order  parameter,  assuming  a 
ferromagnet  for  simplicity.  In  simulated  annealing,  the 
system  is  initialized  at  infinite  temperature,  or  equal  occu¬ 
pation  of  all  classical  spin  states,  and  then  the  temperature  is 
gradually  lowered  to  zero.  The  simulated  spin  dynamics  is 
chosen  to  satisfy  the  detailed  balance  condition  such  that  it 
samples  the  thermal  distribution  at  a  given  (instant)  temper¬ 
ature,  The  initial  state  is  a  paramagnet,  and  therefore  the 
solution,  the  ground-state  spin  configuration  at  zero  temper¬ 
ature,  is  expected  to  have  high  statistical  weight  only  at  low 
enough  temperatures  below  the  ferromagnetic  phase  tran¬ 
sition.  The  exponential  degeneracy  of  the  metastable  state 
corresponds  to  the  entropy  linear  in  the  system  size  N,  which 
significantly  lowers  the  transition  temperature  (see  Fig.  1). 


u 


FIG.  1.  Thermodynamic  functions  of  a  classical  model  char¬ 
acterized  by  the  order  parameter  M  above  (blue  lines)  and 
below  (red  lines)  the  phase  transition  temperature.  Dashed  lines 
correspond  to  potential  energy,  and  solid  lines  correspond  to 
free  energy  that  includes  the  entropy,  shown  separately  by  a 
dash-dotted  black  line.  The  difference  between  dashed  and  solid 
lines  shows  the  contribution  from  entropy.  5 F  is  the  free  energy 
difference  that  has  to  be  overcome  by  classical  Monte  Carlo 
dynamics  at  temperatures  below  the  transition  T  <  Tc. 
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This  can  be  understood  intuitively  from  the  following 
argument.  We  assume  a  mean-field  case  in  which  energies 
of  the  metastable  and  ground  states,  as  well  as  the  barrier 
separating  them,  scale  linearly  with  the  system  size  NEM s, 
NEGS  and  NIJ,  respectively.  We  find  from  equating  the 
free  energies  NQMS  —  NEMS/T  «  —NEGS/T,  where  the 
entropy  imbalance  is  given  by  NQMS,  that  the  transition 
occurs  at  Tc  ~  [(£M§  —  £'gs)/(2ms)]  ~  <5(1)-  Furthermore, 
high  statistical  weight  of  the  ground  state  is  achieved 
only  after  equilibration  at  the  low  temperature  below  the 
phase  transition  T  <  Tc  ~  0(  1).  In  the  presence  of  the 
extensive  barrier  NU,  the  relaxation  towards  the  thermal 
distribution  described  by  the  classical  Kramers  escape 
rate  ~exp(— NU/TC)  is  exponentially  slow.  Moreover, 
the  entropy  gradient  along  the  over-the-barrier  escape  tra¬ 
jectory  gives  rise  to  an  additional  entropic  factor 
~  exp[iV(<2MS  —  Qt)\,  where  (Qr  is  the  entropy  correspond¬ 
ing  to  the  maximum  of  the  free  energy.  Here,  QT  appears  as 
an  exponential  contribution  to  the  prefactor  of  the  Kramers 
rate  [26].  In  this  paper,  we  call  QT  the  “entropic  barrier.” 
The  SA  computation  time  allowing  for  such  relaxation 
to  occur  is  at  least  as  long  as  the  relaxation  time 
tu  ~  exp[yVt/<2MS/ (Ems  —  ■E'gs)  +  N(QMS  —  Qt)\-  1°  fact> 
a  more  careful  analysis  (see  Appendix  A)  shows  that  the 
optimal  SA  computation  time  in  the  presence  of  the  extensive 
entropy  imbalance  is  given  by  the  smallest  of  either  r„  or  the 
exhaustive  search  time  res  ~  2N .  At  the  same  time,  the 
quantum-tunneling  amplitude  saturates  as  T  — ►  0  and  may 
be  more  efficient  than  over  the  barrier  escape,  suggesting  that 
quantum  annealing  could  be  more  efficient  than  simulated 
annealing.  Moreover,  we  see  below  that  the  transverse  field 
introduced  in  QA  lifts  the  degeneracy  of  the  metastable  state, 
and  in  this  way,  QA  avoids  the  entropic  barrier  completely. 
Note,  however,  that  the  quantum-tunneling  rate  in  this  mean- 
field  model  also  scales  exponentially  with  N.  Therefore,  the 
performance  (computation  time)  of  SA  and  QA  is  charac¬ 
terized  by  the  numerical  scaling  factors  in  the  exponent  (in 
front  of  N),  which  have  to  becarefully  compared.  The  result  is 
not  obvious  a  priori  since  here  we  are  comparing  different 
microscopic  mechanisms:  the  quantum  dynamics  con¬ 
strained  by  conservation  laws  with  the  classical  thermal 
excitation  process  constrained  by  the  entropy  imbalance  and 
the  low  temperature. 

The  remainder  of  the  paper  is  organized  as  follows.  In 
Sec.  II,  we  introduce  the  p-spin  model  and  its  WKB  analysis 
describing  the  evolution  of  the  potential  energy  and  the  tran¬ 
sition  rates  in  the  course  of  the  quantum  annealing  algorithm. 
In  Sec.  Ill,  we  discuss  the  dynamics  of  the  model  in  the  course 
of  the  quantum  annealing  and  identify  the  freezing  point  and 
its  optimal  position  in  the  course  of  the  algorithm.  We 
conclude  with  a  discussion  of  the  results  in  Sec.  IV. 

II.  THE  MODEL 

We  consider  N  Ising  spins- 1/2  on  a  fully  connected 
graph  (i.e.,  each  spin  interacts  with  each  other  spin)  with 


uniform  strength  of  spin-spin  and  multispin  interactions, 
subject  to  a  uniform  field.  The  uniform  interactions  give 
rise  to  a  highly  symmetric  Hilbert  space  such  that  the 
unitary  dynamics  of  the  system  is  fully  described  in  terms 
of  the  total  spin  projection  operators  Sa  =  Yl’iLi  of>  where 
a  =  x,  y,  z  and  b°  is  a  set  of  spin- 1/2  operators.  In  other 
words,  the  Hamiltonian  is  defined  as 

H  =  sNf(^Sz\-(l-s)S*.  (1) 

Here,  the  second  term  describes  the  uniform  transverse 
field,  and  the  first  term  is  a  potential  energy  function  f(x), 
which  is  assumed  to  be  a  function  of  the  z-projection 
operator  only.  In  this  paper,  we  consider  a  polynomial  form 
of  f(x),  although  the  general  case  can  be  treated  in  a  similar 
fashion.  Polynomial  terms  of  the  form  /( x)  —  xp  have  a 
natural  microscopic  form  of  p-spin  interaction  of  unit 
strength  H  =  (2 /N)p  of, <3/ ...a]  .  For  example,  a  uni¬ 

form  version  of  the  Sherrington-Kirkpatrick  model 
=  J(I2i°Zi)2  corresponds  to  fix)  =  x2.  In  the 
same  way,  a  fully  connected  graph  with  three-spin  inter¬ 
actions  =  y(JT<7?)3  corresponds  to  the  cubic 

form  of  the  potential  energy  f{x)  =  x3.  Quantum  dynamics 
of  the  fully  connected  graphs  is  not  a  purely  theoretical 
pursuit,  it  can  be  implemented  experimentally  on  existing 
analog  quantum  annealing  hardware  [27],  reducing 
the  model  to  local  interactions  only  [28,29].  Moreover, 
three-spin  interactions  may  be  implemented  directly  in 
superconducting  circuits  [30]  motivated  by  the  novel 
physics  they  could  introduce  [31]. 

Without  loss  of  generality,  we  choose  both  of  the  terms 
in  the  Hamiltonian  Eq.  (1)  to  scale  linearly  with  N.  The 
parameter  s  in  Eq.  (1)  controls  the  relative  strength  of  the 
potential  energy  and  the  transverse  field  and  changes  from 
s  =  0  to  s  =  1  in  the  course  of  the  quantum  annealing 
algorithm. 

For  a  general  function  f(x),  the  Hamiltonian  Eq.  (1) 
commutes  with  S2  =  {Sx)2  +  (Sy)2  +  {Sz)2,  which  is  there¬ 
fore  a  conserved  quantity.  In  the  basis  of  states,  | S,M): 
*S2|S,M)  =  S(S  +  l)|5,Af),  with  definite  total  spin  S  and 
its  projection  on  the  z  axis  M  =  {—.S', . . . ,  .S'} ,  the  matrix 
elements  of  the  Hamiltonian  Eq.  (1)  are  given  by  the 
standard  spin-5  rules, 

SZ\S,M)  =M\S,M),  (2) 

S*f S,M)  =  y/S(S+  1)  —  M[M  ±  1)|5,M  ±  1),  (3) 

where  we  introduced  raising  and  lowering  operators, 
.S'+  (Sx  ±  iSy).  We  introduce  an  integer  parameter 
K  =  0, 1, ....  L7V/2J  to  label  the  total  spin  eigenstates 
S  =  (N /2)  —  K.  The  Hamiltonian  Eq.  (1)  is  symmetric 
with  respect  to  exchanges  of  pairs  of  spins  b,  -o-  bj  and,  in 
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fact,  with  respect  to  all  permutations  of  spins  since  such 
operations  do  not  change  the  sum  over  all  spins  Sa.  This 
symmetry  introduces  high  degeneracy  of  eigenstates  depen¬ 
dent  on  the  total  spin  S.  The  subspace  with  the  maximal  total 
spin  S  =  (N/2)  or  K  =  0  contains  25+1  nondegenerate 
states  (there  are  no  nontrivial  permutations)  corresponding  to 
all  possible  projections  of  the  total  spin  on  the  z  axis.  The 
states  with  K  ±  0  are  highly  degenerate,  with  the  degeneracy 
being  determined  by  the  representations  of  the  group  of 
permutations.  The  eigenstate  with  a  total  spin  labeled  by  K 


has  the  degeneracy 


N 

K  -  1 


~  exp  (NQk), 


where  k  =  (K/N)  =  {0,  (1  /N),  (2 /N), ...,  {l/N)[N/2\}, 
which  corresponds  to  the  entropy  term, 


Qk 


—k\nk  —  (1  —  /r)  ln(  1  —  k)  +  O 


that  has  to  be  added  to  the  free  energy  of  a  state  with  a  given 
energy  E  and  total  spin  parameter  k.T  =  /IE  —  Qk.  Each 
value  of  k  labels  a  Hilbert  subspace  completely  disconnected 
from  that  labeled  by  different  values  of  k.  This  is  a 
manifestation  of  the  spin  permutation  symmetry  that  will 
be  violated  by  the  coupling  to  a  thermal  bath,  which 
introduces  matrix  elements  between  states  with  different  k 
and  effective  relaxation  in  the  system. 

In  this  paper,  we  are  interested  in  a  cubic  potential 
energy, 


rf  \  /  \2  (  ~^max  dmin  \  /  a\ 

f{q)  =  -c{q-q,njn)[q - 2 - )’  t4' 

where  q=(2/N)M,  q  =  {— (1  —  2k), ...,  (1  —  2k)}. 
Equation  (4)  is  the  most  general  cubic  function  with  the 
metastable  minimum  at  qmm  and  the  potential  barrier  top  at 
c/max,  where  f(qmin)  >  /( 1)  ensures  that  q  =  I  is  the  global 
minimum.  Without  loss  of  generality,  we  can  put  c  —  1;  the 
only  effect  of  c  ^  1  is  to  rescale  the  parameter  .s'  in  Eq.  (1). 
We  are  interested  in  analyzing  a  computational  task,  which 
can  be  formulated  for  the  model  Eq.  (1)  by  defining  an 
appropriate  “oracle”  [32]. 

The  cubic  potential  is  chosen  such  that  the  model 
demonstrates  a  first-order  phase  transition  that  is  unavoid¬ 
able  in  the  course  of  quantum  annealing  [34].  Below,  we 
focus  on  the  open-system  quantum  annealing  in  the 
presence  of  dissipation  and  nonzero  temperature,  which 
is  the  case  more  suitable  for  implementation  on  current 
analog  quantum  annealers. 


finding  any  state  within  the  ground-state  potential  well  is 
sufficient  to  find  the  solution  (as  a  local  search  could  allow 
one  to  identify  the  lowest  energy  state  within  the  well). 
Given  the  probability  Vqs  of  finding  the  ground  state  after  a 
single  run  of  duration  v~l,  the  number  of  runs  needed  to 
achieve  this  goal  is  Vq1s.  The  total  computation  time  is 
therefore  given  by 

T  ~  V~l  X  Vqs-  (5) 

The  goal  of  this  paper  is  to  analyze  the  scaling  of  the 
optimal  t  with  the  number  of  qubits  in  the  system  N 
characterized  specifically  by  the  quantity 

Z  =  ^ogT.  (6) 

Below,  we  show  that  the  scaling  of  c  in  our  model  is 
dominated  by  the  exponentially  slow  quantum  tunneling 
through  (or  classical  escape  over)  a  wide  barrier  separating 
the  metastable  and  the  ground  states.  The  tunneling  rate  can 
be  analyzed  using  the  WKB  wave  functions  we  introduce 
in  the  following  sections. 

B.  WKB  wave  functions 

The  dynamics  of  the  model  Eq.  (1)  can  be  described  using 
a  systematic  quasiclassical  WKB  expansion.  For  an  excellent 
review,  see  Ref.  [18].  In  a  spin  model,  this  expansion  is 
performed  in  terms  of  the  small  parameter  e  =  2/N  1, 

which  is  an  analog  of  h  in  the  textbook  WKB  approach  [35]. 
We  consider  a  wave  function  in  the  form 


and  expand  it,  <f>(g)  ~  4>o  +  eTq  +  e2^2  +  G(e2)-  We  can 
further  expand  the  coefficients  4>,  for  a  small  shift  of  the 
argument  from  q  to  q  ±  e, 

^m± i  =  (l  +  o  ±  if® i 

where  O  =  ( dO/dq ).  Substituting  this  expansion  into  the 
Schroedinger  equation,  we  obtain 

where  in  the  main  order  in  e,  the  Hamiltonian  is  diagonal  and 
reads 

e{q)  K  sf(q)  -  .  -(1  -  s)yj (l  -  2k)2  -  q2  cos  4>0.  (7) 


A.  Quantum  annealing  computation  time 

In  the  course  of  the  QA  algorithm,  the  transverse-field 
parameter  .v  is  varied  from  .v  =  0  to  .v  =  1  with  a  fixed  rate  v 
at  a  fixed  inverse  temperature  ji.  The  goal  is  to  find  the 
ground  state  with  a  probability  approaching  1 ,  allowing  for 
repeated  runs  of  the  algorithm.  Here,  we  assume  that 


Note  that  the  function  p  =  4>0  is  precisely  the  canonical 
momentum  conjugate  to  the  coordinate  q,  p  -*  —i(d/ dq).  In 
other  words,  the  second  term  above  is  the  quantum  kinetic 
energy,  which  for  p  <3C  1  corresponds  to  a  particle  with  a 
position-dependent  mass  m~l  =  ( ( I  —  s) \J(  \  —  2kf  —  q2 
moving  in  an  effective  classical  potential, 
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U(k,  q)  =  sf(q)  -  *  (1  -  s)^J (1  -  2k)2  -  q2.  (8) 


Note  that  the  mass  of  the  effective  quantum  particle  increases 
with  increasing  q  and  diverges  as  q  — >  1,  which  affects  the 
efficiency  of  quantum  tunneling  into  the  states  near  q  =  1. 
The  effective  potential  Eq.  (8)  for  different  values  of  s  is 
shown  in  Fig.  2.  In  the  course  of  the  QA  algorithm 
s  —  0  -*  s  —  1,  the  effective  potential  deforms  from  a 
square-root  parabola  corresponding  to  the  ground  state 
(q,k)  =  (0,0)  with  the  maximal  spin  polarization  along 
the  transverse-field  direction  at  5  -►  0  (see  left  panel  in  Fig.  2) 
and  the  classical  potential  corresponding  to  the  ground  state 
( q ,  k)  =  (1,0)  fully  polarized  along  the  axis  of  quantization 
s  — »  1 .  Note  that  the  initial  and  final  states  of  the  algorithm 
are  characterized  by  an  exponentially  small  overlap  (see 
Appendix  B). 

States  with  small  kinetic  energy  are  confined  to  one  of 
the  potential  wells,  centered  around  the  two  minima 

?min(,s)  <  ^lnin(^)  °f  the  effective  potential  U(k,  q),  which 
in  the  course  of  the  evolution  as  s  ->  1,  approach  the 

metastable  q'^n  qmin  and  ground  state  q^n  I  of  the 
classical  model,  respectively.  Confinement  of  the  states  is 
determined  by  the  condition  of  vanishing  classical  velocity, 
v(k,q,E)  =  (dH/dp)  =  0,  which  gives 


-1  <  2 


sf(q)  ~  E 


(l-s)y/(l-2k)2-q2 


<  1. 


(9) 


The  solution  of  this  equation  gives  the  location  of  the 
turning  points  q^p  (k,  E)  and  q-fp  (k,  E),  which  limit  the 


classically  allowed  region.  Inverting  the  secular  equation, 
Eq.  (7),  for  a  given  energy  E,  we  write  the  wave  function  in 
the  leading  order  in  e, 

~  exp  dqp)  >  (10) 

with 


p  =  arccos  2  - 


sf(q)  ~  E 


\l-s)y/(l-2  k)2-qr 
when  the  condition  Eq.  (9)  is  satisfied,  and 

sf(q)  -  E 


p  =  z'arccosh2 


(l-S)^(l-2k)2-q2 


(ID 


(12) 


otherwise.  The  latter  expression  corresponds  to  the  expo¬ 
nentially  decaying  tail  of  the  wave  function  extending 
beyond  the  classically  allowed  region  into  the  potential 
barrier. 


C.  Quantum  phase  transition 

In  the  course  of  the  evolution,  s:  0  -*  1 .  there  is  a  point  of 
zero-temperature  discontinuous  quantum  phase  transition, 
s  =  .VqP|  ,  at  which  the  minimal  energies  in  the  two  wells, 
left  El  at  q  ~  c/min  and  right  ER  at  q  ~  1,  are  equal  to  each 
other.  In  the  course  of  the  QA  algorithm  at  finite  temper¬ 
ature,  this  transition  occurs  at  a  weaker  transverse-field 
strength  s  =  sVY(P)  >  sqPT.  The  phase  transition  point 
can  be  found  from  the  condition  of  equal  occupation  of 


FIG.  2.  The  effective  classical  potential  is  shown  for  5  =  0,  the  pure  transverse-field  case  (left);  s  =  5qPT  «  0.698,  the  zero- 
temperature  quantum  phase  transition  point  (center);  and  5=1,  the  pure  classical  potential  (right).  Different  lines  (bottom  to  top  on  all 

plots)  correspond  to  different  values  of  k  =  0 . 0.5  with  equal  intervals.  Note  that  in  the  absence  of  the  transverse  field,  the  states  with 

different  values  of  k  are  degenerate.  Here,  qmin  =  0  and  qmax  «  0.467. 
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the  two  potential  wells  VL  =VR,  including  the  entropy  of 
the  states.  In  the  large  N  limit  and  at  low  temperatures 
fi  ~  0(N°),  we  can  approximate  the  occupation  number 
'Pl  =  J2E.kexv(-NP-yogZ)Kexp{-NTL-logZ)  by 
a  single  dominant  term  corresponding  to  a  minimum  of 
the  free  energy.  We  write  the  local  minimum  condition  in  a 
potential  well  as  {dJ- /dk)  =  0  to  obtain 


~P 


dE 

Ok 


dQk 

dk 


ln- 


1  -k' 


(13) 


This  determines  the  kmm  corresponding  to  the  minimum  of 
the  free  energy,  whereas  the  optimal  energy  (principal 
quantum  number)  with  a  given  k  corresponds  to  the 
minimum  of  the  potential  U(kmin,q).  We  neglect  the 
quantization  of  levels  due  to  finite  N  for  the  putpose  of 
this  calculation.  At  s  >  Spy(fi).  the  metastable  (left)  poten¬ 
tial  well  is  separated  from  the  ground  state  by  a  potential 
barrier  with  the  shape  determined  by  the  parameters  c/mm 
and  qmax,  and  the  overall  scaling  ~N.  At  low  temperatures, 
the  relaxation  in  this  model  is  therefore  determined  by 
the  rate  of  transitions  between  the  two  wells,  which  we 
calculate  in  the  following. 


with  growing  s  >  s pj(P),  this  state  becomes  metastable. 
The  average  transition  rate  w  W  across  the  barrier  in  the 
presence  of  the  fast  intrawell  relaxation  W  can  be  obtained 
by  calculating  the  total  current  escaping  the  metastable 
well  [38,39], 


WOc-^w(k,E)e  N^E 


(14) 


Equation  (14)  is  a  thermal  average,  weighted  with  the  usual 
Boltzmann  factor  e~NdiE~®k\  including  the  entropy  Qk,  of 
w(k,  E)  ~  | A (k,  E)\2  ~  e~NSv/KB (A~,£) ,  the  incoherent  tunnel¬ 
ing  amplitude  through  the  barrier  of  a  state  with  a  given 
energy  E  and  total  spin  parameter  k.  The  so-called  reduced 
action  5Wkb(^-  E)  can  be  obtained  by  matching  the 
quasiclassical  wave  functions  across  the  banner  region  or 
following  the  analytical  continuation  procedure  [35], 


9* 


SWKB{k,E)  =  -j  dqp(q) 


’<■ 1L 

f  <1r 


2  (sf(q)-E) 


,  dq  arccosh  -  , _ 

Ul  (1  -  s)v/(l  -  2k)2  -  q2 


D.  Quantum  tunneling 

In  a  closed  quantum  system  in  the  absence  of  the  thermal 
bath,  quantum-mechanical  states  are  coherent  superposi¬ 
tion  of  states  in  the  two  wells,  the  so-called  Schroedinger 
cat  states.  Formally,  these  states  correspond  to  a  coherent 
(infinite)  sum  of  amplitudes  of  multiple  tunneling  events 
between  the  wells.  In  a  large  system,  however,  the  level 
splitting  A  corresponding  to  such  a  superposition  state  is 
exponentially  small  (in  the  system  size  N ),  and  therefore, 
such  coherent  dynamics  is  quickly  suppressed  by  small 
perturbations,  such  as  the  hardware  noise.  This  results  in 
overdamped  dynamics  characterized  by  fast  intrawell 
relaxation  towards  thermal  occupation  [36]  reflected  in 
the  level  width  W  »  A  and  exponentially  rare  incoherent 
tunneling  events  with  the  rate  ~A2.  This  is  the  regime  we 
consider  in  this  paper.  At  the  same  time,  we  neglect  the 
effect  of  noise  on  the  tunneling  event  itself  since  tunneling 
is  a  fast  process  occurring  on  the  time  scale  1  / Q,  where  Q 
is  the  frequency  determined  by  the  curvature  of  the 
potential.  We  assume,  therefore,  that  £2  W  A.  We 
are  interested  only  in  the  exponential  scaling  of  the 
transition  rates  in  this  paper,  ignoring  the  renormalization 
of  the  preexponential  factor  that  may  be  substantial  in  the 
regime  of  strong  coupling  to  the  environment.  Note  that 
the  overdamped  dynamics  and  thermalization  are  expected 
even  in  the  absence  of  the  coupling  to  a  thermal  bath;  it  can 
be  introduced  by  a  weak  disorder  in  the  spin-spin  inter¬ 
actions  8H  =  e{jd]aZj  or  even  a  weak  random  transverse 
field  [37], 

At  s  <  Spj(fi),  the  ground  state  corresponds  to 
q  =  <7min  (>v) -  As  the  system  goes  past  the  phase  transition 


The  sum  in  Eq.  (14)  can  be  approximated  by  its  largest 
term,  the  rest  being  exponentially  smaller, 

w  ~  —  ^\v(k,  E)e~N^E~ ~  e~NS,  (15) 

Z  k,E 


where  the  largest  term  is  found  by  minimizing  the  action, 
S  =  SWKB+0E-Qk  +  ±lnZ  (16) 

with  respect  to  k  and  E, 


r(£)  =  _(«wp  =  ^  (17) 


^wkb  qdE  _  ()Qk  ^  k 
dk  dk  dk  1  —  k 


(18) 


In  Eq.  (17),  fixed  k  is  assumed.  Since  Qk  is  independent  of 
the  energy  level  E,  the  conditions  on  the  optimal  tunneling 
parameters  separate  into  the  standard  condition  Eq.  (17) 
[requiring  the  period  of  motion  in  the  inverted  potential 
T(E)  to  match  the  inverse  temperature  [i  [39]]  and  the 
condition  Eq.  (18)  due  to  the  entropy  of  states  dependent  on 
k,  which  introduces  novel  physics  in  the  dynamics  of  this 
model.  Equations  (16)— (18)  need  to  be  supplemented  with 
conditions  ensuring  that  the  energy  and  total  spin  k  are 
conserved  in  the  tunneling  event  (we  emphasize  that  we 
neglect  here  the  effect  of  the  thermal  bath  during  the 
tunneling  event).  The  tunneling  from  the  metastable  well 
has  to  be  at  energy  and  spin  values,  E  and  k,  at  which  a  state 


021028-6 


OPEN-SYSTEM  QUANTUM  ANNEALING  IN  MEAN-FIELD  ... 


PHYS.  REV.  X  6,  021028  (2016) 


FIG.  3.  Tunneling  action  as  a  function  of  inverse  temperature  /?, 
for  q ^  =  0,  qmax  =  0.533  at  s  =  0.85  >  Sqpt-  Left:  Black  dots 
and  the  solid  blue  line  fit  correspond  to  an  optimal  finite- 
temperature  quantum-tunneling  action.  The  horizontal  dotted 
line  is  the  p  — ►  oo  limit  of  the  incoherent  tunneling  action 
corresponding  to  the  quantum-mechanical  tunneling  from  the 
lowest  level  of  the  metastable  well  corresponding  to  k  =  0.  The 
solid  black  line  corresponds  to  the  optimal  action  of  the  classical 
Glauber  dynamics.  Blue  and  red  dashed  lines  show  the  classical 
over-the-barrier  escape  action  at  k  =  0  and  along  the  line  q  = 
1  —  2k  with  the  potential  maximum  at  kt  «  0.152,  respectively. 
The  latter,  kt,  is  the  point  of  inflection,  i.e.,  the  point  where  the 
minimum  of  the  right  potential  well  merges  with  the  maximum  at 
the  barrier  top  (see  Fig.  2).  The  vertical  dashed  line  corresponds 
to  the  temperature-driven  phase  transition  point  pPT « 4.32 
corresponding  to  equal  occupation  of  the  two  potential  wells. 
Inset:  T(E)  =  —dS/dE  at  fixed  k,  which  corresponds  to  the 
period  of  quantum-mechanical  tunneling  trajectory  in  the  imagi¬ 
nary  time  representation.  The  dominant  contribution  comes  from 
tunneling  at  the  energy  determined  from  T(E)  =  P  for  p  >  n; 
at  P  <  Tmin,  the  transition  rate  is  dominated  by  the  over-the- 
barrier  escape.  Different  lines  from  bottom  to  top  correspond  to 
different  values  of  j  at  fixed  k  =  0.  Red  dots  correspond  to  Tmm. 
Note  that  the  minimum  Em in:  T(Emin)  =  Tmin  does  not  corre¬ 
spond  to  the  highest  tunneling  energy.  This  means  that  the 
transition  from  the  quantum-tunneling  regime  to  over-the-barrier 
escape  has  discontinuous  first-order  character.  Right:  Quantum- 
classical  transition  region  on  a  larger  scale.  Black  dots  correspond 
to  the  dominant  tunneling  action  at  each  p.  Different  lines  show 
the  quantum-mechanical  tunneling  action  for  fixed  k\  colors 
correspond  to  growing  0  <  k  <  0. 152  (red  to  blue).  The  quantum¬ 
tunneling  process  conserves  energy  and  the  total  spin  value  k. 
At  k  >  kLR,  the  point  where  EL(kLR)  =  ER(kLR),  the  quantum 
mechanical  tunneling  process  requires  the  state  in  the  metastable 
well  to  be  at  an  energy  E>ER{k),  which  comes  with  an 
additional  thermal  excitation  cost  P[ER (k)  —  EL(k)\  in  the  tun¬ 
neling  action.  Therefore,  for  growing  k,  the  tunneling  action  S{p) 
resembles  linear  classical  dependence.  Solid  black  and  dashed 
red  and  blue  lines  are  the  same  as  in  the  left  figure. 

exists  in  the  ground  state  well,  i.e.,  E  >  min  {ER(k)}, 
which  is  not  always  satisfied  in  the  system  with  large 
entropy  of  states;  i.e.,  T L  <  TR  does  not  necessarily 
imply  El  <  Er. 

Equation  (17)  has  a  solution  in  a  range  of  energies  E 
such  that  Tmin  <  T(E)  <  oo  (see  inset  in  Fig.  3,  left).  In  the 
case  of  [k  >  Tmin.  the  quantum-tunneling  process  dominates 
in  the  sum  in  Eq.  (15).  For  p  <  Tmm,  there  are  no  solutions 


to  Eq.  (17),  and  therefore,  the  optimal  energy  is  at  the  edge 
of  the  interval  E  =  U (qmax)  —  U (qmm)  corresponding  to 
the  height  of  the  barrier.  In  other  words,  in  this  regime,  the 
over-the-barrier  escape  process  dominates,  with  p  ~  Tmm 
being  the  point  of  a  quantum-to-classical  phase  transition. 
Note  that  the  global  minimum  Emm  of  the  function  T (E)  in 
the  inset  in  Fig.  3  does  not  always  correspond  to  the  top  of 
the  barrier,  which  means  that  the  quantum-to-classical 
transition  (in  the  limit  N  -*■  oo)  has  a  discontinuous  first- 
order  character  [19].  Considering  Eq.  (18),  we  look  for  a 
solution  in  the  interval  0  <  k  <  kt,  where  k,  is  the 
inflection  point  of  the  potential  U(q.k)  in  which  the  right 
(the  ground-state)  potential  well  disappears  (see  Fig.  2); 
since  quantum  tunneling  conserves  the  total  spin  k,  for  it  to 
occur  there  must  exist  states  with  matching  k  in  the  ground- 
state  potential  well. 

The  result  of  this  optimization  procedure  is  the  optimal 
action  S(JJ)  at  a  fixed  s  shown  in  Fig.  3.  In  the  vanishing 
temperature  limit  p  oo.  where  the  effect  of  entropy  on 
the  occupation  of  levels  is  negligible,  S(p)  corresponds  to 
the  quantum  tunneling  from  the  lowest  energy  level  in  the 
metastable  well  corresponding  to  k  =  0  (horizontal  dashed 
line  in  Fig.  3,  left).  As  temperature  increases  (/)  decreases), 
the  entropy  starts  playing  a  role  in  dynamics,  and  S(P) 
increases  up  to  some  maximum  value  (see  blue  solid  line  in 
Fig.  3,  left),  in  contrast  to  the  usual  (nondegenerate)  case 
[39]  where  the  quantum-tunneling  rate  increases  monoto¬ 
nously  with  increasing  temperature  /T  1 .  This  is  a  result  of 
the  entropy  providing  high  statistical  weight  to  the  state 
with  a  suboptimal  tunneling  rate.  This  behavior  is  natural 
because  the  transverse  field  splits  the  degeneracy  of  the 
metastable  state  at  q  =  q[mn  favoring  the  state  with  maxi¬ 
mal  total  spin,  corresponding  to  k  =  0.  At  the  same  time, 
the  transverse  field  provides  maximal  quantum  kinetic 
energy  (minimal  mass  in  the  quadratic  approximation)  to 
the  same  state,  as  can  be  seen  from  Eq.  (7),  whereas  the 
finite  temperature  due  to  the  effect  of  entropy  favors  the 
states  with  1/2  >  k  >  0  (at  p  ->  0),  which  have  lower 
kinetic  energy  (higher  mass)  and  correspond  to  the  sub- 
optimal  tunneling  rate  through  the  classical  potential  sf(q). 
The  result  is  S(P)  increasing  with  decreasing  p  up  to  the 
regime  of  the  transition  into  the  classical  escape  at 
high  temperatures  p  <s:  7j,lln.  The  classical  over-the-barrier 
excitation  is  described  by  the  thermal  excitation  rate 
wcl  ~  Z  1  exp [—N(pE  —  Qp)].  Note  that  the  classical  proc¬ 
ess  is  driven  by  the  Glauber  dynamics  of  the  spins  due 
to  the  effect  of  the  thermal  bath.  This  process  does  not 
conserve  the  total  spin  k,  and  therefore,  the  optimal 
classical  trajectory  is  determined  by  the  saddle  point  of 
the  free  energy  (including  the  entropy)  in  the  2D  space  of 
( k,q ).  Entropy  provides  an  additional  cost,  reducing  the 
transition  rate,  which  is  reflected  in  the  finite  offset  of  the 
dependence  of  the  classical  transition  rate  on  p  as  p  -*  0 
(see  black  solid  line  in  Fig.  3,  left;  see  Appendix  A  for  more 
details).  Blue  and  red  dashed  lines  in  Fig.  3  indicate  the 
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linear  in  ji  dependence  of  the  energy  cost  of  the  over-the- 
barrier  excitation  at  k  =  0  and  along  the  line  of  q  —  1—2 k 
with  the  potential  maximum  at  £  «  [the  latter  is  the 
inflection  point  of  the  potential  U(k,q)],  respectively, 
which  correspond  to  the  low-temperature  potential- 
energy-dominated  regime  and  the  high-temperature 
entropy-dominated  regime  [40]. 

The  steepest  descent  approximation  Eq.  (15)  is  appli¬ 
cable  as  long  as  the  preexponential  factors  in  the  sum  are 
nondivergent,  which  is  true  away  from  ji  ~  Tmm,  a  phase 
transition  point  from  the  quantum-tunneling  regime  to  the 
classical  over-the -barrier  escape  regime  [39].  On  general 
grounds,  we  expect  the  action  to  be  continuous  even  in  the 
case  where  this  quantum-to-classical  transition  is  of  the  first 
order,  with  the  discontinuity  occurring  in  the  derivative  of 
the  action  [19].  Therefore,  we  expect  Fig.  3  to  provide  a 
qualitatively  correct  dependence  of  S(ft)  in  the  whole  range 
of  inverse  temperatures  j>. 

III.  QUANTUM  ANNEALING 
COMPUTATION  TIME 

We  now  turn  to  the  calculation  of  the  computation  time 
defined  in  Eq.  (6).  The  quantum  mechanical  tunneling  rate 
vanishes  as  s  -*  1  (at  very  low  temperatures  of  interest, 
here  the  over-the-barrier  transition  rate  is  also  weak). 
Therefore,  there  exists  a  point  s  =  sF  in  the  course  of 
the  sweep  of  the  transverse  field  where  the  relaxation  time 
~w— 1  (s/;)  required  to  achieve  thermal  distribution  becomes 
longer  than  the  length  of  the  algorithm,  vc~ 1  (sF)  ~  v~l.  In 
other  words,  the  computation  will  be  finished  before  the 
thermal  equilibrium  is  reached.  The  system  effectively 
freezes  the  values  of  the  occupation  numbers  of  the  left 
VL  and  right  VR  wells  at  the  last  point  (in  the  course  of  the 
algorithm)  where  the  interwell  transition  process  was  still 
fast  enough  (the  intrawell  relaxation  may  still  be  efficient  at 
s  >  sF).  We  call  this  the  freezing  point.  Sudden  freezing  of 
quantum  dynamics  in  the  course  of  quantum  annealing  is 
typical  for  the  algorithms  relying  on  quantum  tunneling 
between  states  separated  by  large  Hamming  distance 
[3,41-44].  The  computation  time  of  the  algorithm  (its 
exponential  scaling)  is  therefore  determined  by  the  thermal 
(equilibrium)  occupation  probability  'PGS  =  VR(sF),  and 
the  equilibration  time  vv~ 1  (sF)  at  the  freezing  point, 

£  =  ^logr  ~  ^l°g  [gJVSU  +  eN^R~^)],  (19) 

where  S  is  the  optimal  quantum-tunneling  action  given 
by  Eq.  (16),  with  s  =  sF.  The  optimal  computation  time 
of  the  quantum  annealing  is  found  by  minimizing  with 
respect  to  the  location  of  sF  and  the  temperature  ji  at  which 
the  computation  is  performed.  The  point  of  the  phase 
transition,  [i  =  /)PX  and  sF  =  .syr,  respectively,  defined 
by  !F R  —  T L  =  0,  separates  two  scaling  regimes  of  the 
computation  time, 


»  $  +  P  <  Pvr(s) 

VS  P>Pw{s). 

The  high-temperature  limit  (Ji  — >  0)  of  this  expression  is 
given  by  the  entropy  difference  between  the  two  wells.  This 
is  the  limit  of  a  local  exhaustive  search,  which  is  a  bound  on 
the  computation  time  of  the  algorithm,  and  we  need  to 
compare  it  to  the  optimal  computation  time  £(/?)  to  make 
sure  we  find  the  global  optimum.  We  first  consider  [i  <  /iPX, 
assuming  the  free  energy  of  the  system  demonstrates  two 
minima,  T R  and  T L.  In  the  quantum-tunneling  regime, 
both  S(Ji)  and  J-  R  —  J-  L  decrease  with  growing  /?,  and 
therefore,  £(/))  is  monotonously  decreasing  as  well.  In  the 
classical  regime,  S  +  (F  R  —  T L  is  also  a  monotonous 
function,  which  depending  on  competition  between  S(jj) 
and  T R  —  J- L,  can  be  either  increasing,  in  which  case  j  — ►  0 
is  the  optimal  classical  computation  regime  (i.e.,  the  local 
exhaustive  search  limit),  or  decreasing  towards  the  critical 
point  [i  =  /iP-|  .  Therefore,  we  need  to  analyze  the  perfor¬ 
mance  in  the  regime  / 3  >  [ivv  and  s  >  sPX  and  compare  it  to 
a  local  exhaustive  search.  The  optimal  computation  time  in 
this  regime  [see  Eq.  (20)]  is  determined  by  the  minimum  of 
the  transition  action  S  with  respect  to  fi  and  s. 

Because  of  the  concave  dependence  of  S(/i)  on  the 
inverse  temperature,  as  shown  in  Fig.  3  (left),  the  minimum 
of  the  action  with  respect  to  [i  corresponds  to  one  of  the 
edges  of  the  inverse  temperature  interval,  i.e.,  either 
P  =  [i\ry  or  fi  — >  co.  The  global  minimum  is  therefore  the 
smallest  of  5|/?PX(sf'),  sF]  and  S(oo,sF).  The  minimum 
with  respect  to  sF  (and  therefore  the  global  minimum)  is 
determined  by  comparing  these  two  functions. 

A  typical  critical  line  is  shown  in  the  left  panel  of  Fig.  4. 
The  inverse  critical  temperature  pPT(s),  blue  solid  line  in 
Fig.  3  (left),  diverges  at  the  point  of  the  quantum  phase 
transition  in  the  course  of  the  algorithm, =  .vqPX  (vertical 
dashed  line  in  Fig.  4),  and  at  s  >  5qPX  it  monotonously 
decreases,  with  s  approaching  the  classical  transition 
temperature  at  s  =  1  (horizontal  dashed  line  in  Fig.  4). 
The  right  panel  of  Fig.  4  shows  the  over-the-barrier  escape 
action  at  the  given  critical  temperature  fipj(s),  blue  (upper) 
solid  line.  The  blue  dots  in  the  left  and  right  panels 
correspond  to  the  same  values  of  s.  The  classical  action, 
following  the  behavior  of  the  inverse  critical  temperature, 
diverges  at  the  point  of  the  quantum  phase  transition  s  = 
SqPX  and  decreases  with  s  in  the  course  of  the  algorithm 
approaching  the  value  corresponding  to  the  classical  model 
s  =  1,  shown  as  the  blue  (upper)  dashed  line.  The  latter 
is  the  optimum  (due  to  the  highest  critical  temperature) 
and  therefore  defines  the  optimal  computation  time  of  the 
classical  simulated  annealing  algorithm  (see  Appendix  A 
for  more  details).  The  optimum  of  the  quantum  action,  at 
vanishing  temperature  /?  -*■  oo,  is  given  by  the  tunneling 
from  the  bottom  of  the  metastable  well.  In  this  limit,  the 
entropy  does  not  affect  the  occupation  of  the  energy  levels 
in  the  course  of  the  algorithm.  We  assume  that  the 
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FIG.  4.  Left  panel:  Inverse  critical  temperature  dependence  on 
the  transverse-field  strength  /tpj/s)  for  gmin  =  0.88,  qmm  =  0.955. 
The  horizontal  dashed  line  corresponds  to  the  classical  transition 
temperature,  and  the  vertical  dashed  line  marks  the  point  of 
the  quantum  phase  transition.  Right  panel:  Comparison  of  the 
vanishing  temperature  quantum-tunneling  action  to  that  of  the 
classical  over-the  barrier  escape.  The  blue  line  shows  optimal 
over-the-barrier  escape  action  along  the  critical  lineshown  in  the 
left  panel.  Blue  points  correspond  to  the  same  values  of  s  as  in  the 
left  figure.  The  horizontal  blue  line  corresponds  to  the  transition 
rate  at  the  critical  temperature  in  the  classical  model,  $  =  1. 
The  vertical  line  corresponds  to  the  point  of  zero-temperature 
quantum  phase  transition.  The  red  line  corresponds  to  the 
vanishing  temperature  limit  of  the  quantum-tunneling  action. 
The  red  horizontal  dashed  line  corresponds  to  the  minimum 
of  the  quantum-tunneling  action  at  the  quantum  phase  transition 
point.  The  local  exhaustive  search  corresponds  to 
|  =  (1  /N)  logres  =  Qci(qaan)  a  0.227,  larger  than  the  quantum 
annealing  time. 

temperature  is  still  high  enough  such  that  there  is  suffi¬ 
ciently  fast  intrawell  relaxation.  Figure  4  (right  panel) 
shows  S(oo,s),  red  (lower)  solid  line,  which  assumes  the 
minimal  value  at  the  quantum  phase  transition  point 
s  =  sqpt-  Note  that  this  is  not  true  for  all  the  parameter 
values;  instead,  the  minimum  of  S(oo,s)  often  occurs  at 
some  ^  >  sqpt,  and  its  value  at  this  point  can  be  as  much  as 
2  times  smaller  than  the  value  at  s  =  SqPX  (see  Fig.  5, 
left  panel). 

Figure  4  (right  panel)  demonstrates  that  the  quantum¬ 
tunneling  process  may  be  more  efficient  than  over-the- 
barrier  escape.  Both  classical  and  quantum  transition  rates, 
and  therefore  the  corresponding  computation  times,  scale 
exponentially  with  the  system  size,  r  oc  exp(cdV),  yet  the 
coefficient  in  the  exponent  is  smaller  in  the  case  of  QA 
as  compared  to  the  classical  simulated  annealing,  which 
corresponds  to  a  polynomial  speedup.  Note  that  it  is 
important  to  compare  the  computation  times  in  Fig.  4  with 
that  of  the  local  exhaustive  search  in  the  interval  (gmin,  1), 
which  is  the  high-temperature  limit  of  SA.  The  correspond¬ 
ing  computation  time  is  given  by  the  entropy  £  =  Qci(Vmin)> 
which  for  parameters  considered  in  Fig.  4,  <2cl(0.88)« 
0.227,  scales  worse  than  the  QA  computation  time  shown 
by  red  dots  in  Fig.  4  (right  panel).  We  further  compare 
the  optimal  performance  of  the  open-system  quantum 
annealing  Sopt  to  simulated  annealing  for  a  wide  range 


FIG.  5.  Left  panel:  Quantum  mechanical  tunneling  action  at 
vanishing  temperature  in  the  incoherent  regime  as  a  function  of 
the  transverse-field  parameter  s.  Different  curves  correspond  to 
different  values  of  gmax  =  0.929,  0.946,  0.958,  0.961  (order  is 
from  bottom  to  top  on  the  right  side  of  the  plot)  with  fixed 
?mm  =  0-9-  Note  that  the  optimal  tunneling  rate  emphasized  by 
the  horizontal  dashed  lines  does  not  always  correspond  to  the 
phase-transition  point  denoted  by  vertical  dashed  lines.  Right 
panel:  Optimal  action  for  SA  (blue)  and  QA  (red)  at  gmin  =  0.88 
for  different  values  of  qmax  [see  Eqs.  (19)  and  (20)  and  the  text 
for  details].  The  dashed  red  line  corresponds  to  the  coherent 
quantum-tunneling  action  (scaling  as  1/2  of  the  action  in  the 
incoherent  regime).  The  vertical  dashed  line  corresponds  to 
?min  =  i  <?max>  at  which  point  the  classical  potential  has  a 
degenerate  ground  state  at  q  =  gmin  and  q  =  1 .  At  this  point, 
both  the  SA  action  and  the  QA  action  diverge;  however,  the  QA 
action  diverges  logarithmically  slow,  in  contrast  to  the  SA  action. 
For  the  parameters  chosen,  the  local  exhaustive  search  corre¬ 
sponds  to  £  =  (1  /N)  log  res  =  Qd  (qm in)  «  0.227,  which  is  above 
the  value  of  the  QA  action  in  the  figure,  except  for  the  close 
vicinity  of  the  divergence  point,  i.e.,  QA  is  faster  than  SA  and  an 
exhaustive  search  for  most  of  the  parameter  values. 

of  barrier  shapes  within  the  cubic  model  Eq.  (4)  by  varying 
the  location  of  the  metastable  minimum  c/min  and  the  barrier 
top  r/milx  in  Eq.  (4)  (see  right  panel  in  Figs.  5  and  6).  QA 
outperforms  S  A  in  a  range  of  the  parameter  space  where  the 
potential  barrier  separating  the  metastable  and  the  ground 
states  is  small.  Note  that  the  origin  of  the  speedup  in  the 
case  considered  here  is  distinct  from  the  standard  intuition 
of  thin  and  tall  barriers  favoring  quantum  tunneling  (see 
Ref.  [13]  for  refining  the  standard  intuition)  since  the  shape 
of  the  potential  is  cubic  throughout  the  range  of  parameters 
shown  in  Fig.  6.  Instead,  the  quantum  algorithm  turns 
out  to  be  more  efficient  because  it  proceeds  along  a  path 
with  lower  entropic  cost  than  the  path  that  SA  takes.  This 
is  a  result  of  the  transverse  field  lifting  the  degeneracy  of 
the  metastable  state  in  the  quantum  case.  Therefore,  the 
smallness  of  the  barrier  required  for  the  speedup  in  this  case 
is  determined  by  the  comparison  of  the  quantum-tunneling 
action  and  the  SA  action,  including  entropy  cost  of  over- 
the-barrier  escape,  the  latter  being  a  combination  of  fipjU 
and  the  additional  entropic  cost  Qci(c/mm)  —  <2c|(<yUip), 
where  gtop  is  the  maximum  of  the  free  energy  (see 
Appendix  A  for  details).  Note  also  that  the  numerical 
value  of  the  ratio  of  the  logarithms  of  the  normalized 
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FIG.  6.  Left  panels:  Potential  barrier  between  the  metastable  and 
the  ground  state,  Eq.  (4),  corresponding  to  s  =  1  for  two  pairs  of 
values  (<7min.  <7max)  equal  to  (0,  0.633)  (upper  plot)  and  (0.88, 
0.955)  (lower  plot).  Right  panel:  Scaling  exponent  of  the  optimal 
computation  time  of  quantum  annealing  (red)  [see  Eqs.  (6)  and 
(19)]  and  simulated  annealing  (blue)  given  by  ( 1  / N )  log  r,  with  z 
given  by  Eq.  (A2)  as  a  function  of  the  location  of  the  metastable 
minimum  c/min  and  the  top  of  the  barrier  gmax.  The  lower  of  the  two 
surfaces  corresponds  to  the  shorter  computation  time.  QA  is 
advantageous  for  a  range  of  parameters  corresponding  to  suffi¬ 
ciently  narrow  potential  barriers.  Solid  lines  in  the  plane  S  =  0 
correspond  to  q^m  =  gmax  and  qmin  =  |  qmm  (under  both  surfaces 
in  the  figure)  outlining  the  range  of  possible  potentials  in  the  cubic 
model  Eq.  (4). 

computation  times  £  =  (l/N)  log  r  for  QA  and  SA  can  be 
substantial  [see  Fig.  5  (right  panel)]. 

IV.  DISCUSSION 

In  this  paper,  we  consider  a  model  problem  of  open- 
system  quantum  annealing  that  allows  analytical  inves¬ 
tigation  and,  at  the  same  time,  demonstrates  key  features 
of  complex  optimization  problems,  including  the  discon¬ 
tinuous  first-order  phase  transition  and  the  exponential 
degeneracy  of  the  metastable  state.  We  demonstrate  that, 
for  problems  with  extensive  degeneracy,  the  SA  algorithm 
relies  on  over-the -barrier  escape  at  very  low  temperatures, 
0(1);  as  a  result,  the  SA  computation  time  is  expo¬ 
nential  in  N.  In  addition,  the  classical  over-the-barrier 
escape  rate  is  reduced  because  of  an  additional  entropic 
barrier,  the  difference  in  the  statistical  weights  of  the 
metastable  well  and  the  states  near  the  top  of  the  barrier. 
At  the  same  time,  we  show  that  a  computational  advantage 
can  be  gained  using  open-system  quantum  annealing, 
which  exploits  the  effects  of  thermally  assisted  quantum 
tunneling  and  relaxation.  In  the  course  of  QA,  the  applied 
transverse  field  splits  the  degeneracy  of  the  metastable 
state  and  provides  the  nondegenerate  low-energy  states 
with  the  highest  quantum  kinetic  energy;  in  other  words, 
the  transverse  field  favors  the  quantum-tunneling  trajecto¬ 
ries  involving  nondegenerate  states.  In  this  way,  the  QA 
algorithm  naturally  circumvents  the  entropic  barrier  and 
thus  gains  a  clear  advantage  over  SA.  This  is  a  novel  feature 
in  the  quantum-tunneling  process  caused  by  the  entropy  of 


the  metastable  states,  manifested  particularly  in  the  tunnel¬ 
ing  rate  decreasing  with  increasing  temperature.  As  a 
consequence,  the  optimal  quantum  annealing  regime  cor¬ 
responds  to  vanishing  temperature;  i.e.,  raising  the  temper¬ 
ature  reduces  the  efficiency  of  QA  [45].  The  faster 
tunneling  with  lower  temperature  contrasts  the  standard 
intuition  that  raising  temperature  should  always  improve 
transition  rates.  In  this  paper,  we  provided  an  important  and 
very  widespread  counterexample  to  this  standard  intuition. 
We  also  found  that  at  low  temperatures,  the  optimal 
quantum-tunneling  rate  does  not  always  correspond  to 
the  point  of  the  phase  transition  SqPX;  in  fact,  tunneling 
at  s  >  SqPT  can  have  a  substantially  higher  rate,  which  can 
be  exploited  in  conjunction  with  the  noise-induced  thermal- 
ization  to  improve  the  performance  of  the  algorithm.  In  the 
mean-field  model  considered  in  this  paper,  the  comparison 
of  the  quantum  annealing  and  simulated  annealing  comes 
down  to  the  numerical  coefficient  in  the  scaling  of  the 
computation  time  with  N.  We  demonstrated  that  optimal 
QA  could  outperform  SA  in  a  certain  parameter  range  of 
our  model  characterized  by  small  potential  barriers.  This  is 
in  spite  of  the  constrained  nature  of  the  quantum-tunneling 
process  due  to  the  conservation  laws,  as  opposed  to  the 
unconstrained  classical  Glauber  dynamics  of  SA. 

We  also  identified  key  features  of  the  model  affecting  the 
performance  of  the  QA,  specifically,  the  quantum  fluctua¬ 
tions  strength  in  the  vicinity  of  the  phase  transition  point  and 
the  diverging  mass  in  the  quantum  kinetic  energy,  which 
both  strongly  affect  the  efficiency  of  quantum  tunneling  in 
the  course  of  the  algorithm.  In  this  respect,  it  would  be 
interesting  to  use  methods  developed  in  this  paper  to  explore 
QA  in  mean-field  models  with  different  types  of  driver 
Hamiltonians  whose  ground  states  are  not  simple  product 
states  where  ferromagnetic  order  competes  with  the  trans¬ 
verse  (XT)  ferromagnetism  or  superfluidity  [ 46 — 49  ]. 

Note  that  the  model  of  escape  from  a  highly  degenerate 
potential  well  considered  in  this  paper  is  quite  generic,  as  it  is 
natural  to  expect  that  the  applied  transverse  field  will  favor 
the  trajectories  involving  nondegenerate  states,  thus  circum¬ 
venting  the  entropic  barrier,  which  slows  down  classical 
stochastic  dynamics.  This  is  because  the  states  most  strongly 
affected  by  the  transverse  field  both  acquire  the  largest 
energy  shifts  and  the  highest  kinetic  energy  and  hence  the 
highest  transition  rates.  This  suggests  a  speculation  that 
optimization  problems  in  which  entropy  is  a  dominant  factor 
yet,  at  the  same  time,  which  are  characterized  by  a  potential 
energy  landscape  that  can  be  exploited  for  a  more  efficient 
search  (Grover’s  unstructured  search  problem  [50]  with  an 
added  potential  landscape),  may  represent  a  class  of  prob¬ 
lems  where  quantum  annealing  could  have  a  computational 
advantage  over  simulated  annealing. 
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APPENDIX  A:  SIMULATED  ANNEALING 

We  discuss  the  performance  of  simulated  annealing  on 
the  model,  Eqs.  (1)  and  (4),  in  the  purely  classical  case  of 
5  =  1  (see  the  rightmost  panel  in  Fig.  2).  The  SA  algorithm 
is  realized  by  starting  with  the  infinite  temperature  limit 

-*  0,  i.e.,  equal  occupation  of  all  states,  and  reducing 
the  temperature  to  zero.  For  simplicity  (and  without  loss 
of  generality),  we  consider  an  algorithm  where  (i  changes 
linearly  in  time  from  0  to  a  very  large  value  with  a  fixed 
rate  v. 

The  time  it  takes  to  find  the  ground  state  with  probability 
approaching  around  1,  allowing  for  repeated  runs  of  the 
algorithm,  is  given,  as  in  the  quantum  case,  by 

Simulated  annealing  relies  on  the  system  reaching  thermal 
equilibrium  throughout  at  least  part  of  the  algorithm,  such 
that  the  ground-state  occupation  is  given  by  the  Gibbs 
distribution.  In  problems  with  a  free  energy  barrier  sepa¬ 
rating  the  initial  and  the  ground  state,  the  system’s 
relaxation  time  is  dominated  by  the  over-the-barrier  escape 
probability  with  T{q)  =  fiE(q)  +  Qc\(q), 

w{fi)  ~  exp  [-N (qtop)  -  F{qm in))], 

given  by  the  statistical  weight  of  the  escape  trajectory  [26], 
where  we  include  the  entropy  of  the  classical  state  with 
magnetization  q  given  by  Qcl(<?)  ~  — [(1  +  r?)/2]  ln[(l-|- 
q)/2]  —  [(1  —  q)/2]  ln[(l  —  q)/2],  and  qtop  is  the  maxi¬ 
mum  of  T (q).  Here,  w(JT)  reduces  exponentially  with 
decreasing  temperature  (growing  ff),  and  therefore,  in 
analogy  with  the  quantum  case  considered  in  the  main 
text,  there  exists  a  freezing  point  [)  =  ()F  in  the  course  of  the 
sweep  of  the  inverse  temperature  where  the  relaxation  time 
~w_1  (/i/,-)  required  to  achieve  thermal  distribution  becomes 
longer  than  the  length  of  the  algorithm,  w~l{pF)  ~  v~l . 
After  this  point,  the  values  of  the  occupation  numbers  of  the 
left  Vj  and  right  V R  potential  wells  are  effectively  frozen. 
Therefore,  the  computation  time  is  determined  by  two 
quantities  calculated  at  the  freezing  point,  the  occupation 
of  the  ground-state  potential  well  'PGS  =  V  R{fiF)  and  the 
relaxation  time  at  the  freezing  point  w~l{f)F), 

T  ~  £'VU^topU‘U<7IIlm)]  (1  —  -^Urnm)]  (Al) 

where  we  keep  only  the  main  order  in  the  limit 
N  ->  oo  such  that  VGs{PF)  ~  exp[— jF(l)]/[exp(— JF(1)]  + 
exp[— J-{qmm)]-  The  optimal  computation  time  can  be 
obtained  by  minimizing  with  respect  to  the  inverse  temper¬ 
ature  at  the  freezing  point,  (d/dpF)[(l/N)\og2r\  =  0.  This 


derivative  is  discontinuous  at  the  point  of  the  phase  transition 
/ipx,  where  J-R  —  T L  —  0.  The  computation  time  is  an 
increasing  function  at  [)F  >  /?PT,  as  it  is  dominated  by  the 
decreasing  transition  rate,  the  prefactor  in  front  of  the  curly 
brackets  in  Eq.  (Al).  However,  it  can  be  either  a  monoto¬ 
nously  decreasing  or  increasing  function  at  fiF  < 
depending  on  the  competition  between  the  prefactor  and 
the  exponents  in  the  brackets  in  Eq.  (Al).  Therefore,  the 
global  minimum  of  r(/)F)  corresponds  to  the  smallest  value 
out  of  r(/?PT)  and  t(0)  «  2N,  in  the  decreasing  and  increasing 
cases,  respectively.  The  latter  corresponds  to  the  exhaustive 
search,  i.e.,  2N  repetitions  of  infinitely  fast  SA.  The 
ground  state  at  q  =  I  is  unique,  Q{\)  ~  0;  therefore, 
the  point  of  the  classical  phase  transition  is  at 
Art  =  [<2(4min)/£(?min)  -£(!)]•  Therefore,  the  optimal 
SA  computation  time  corresponds  to  the  smaller  of  the 
two  values, 


exp 

\7  f  top) 

L  v 

^^Q(qmin)+SQ) 

2N, 

where  5Q  =  Q(qmm)  —  Q{qtov)-  Note  that  the  entropy  asso¬ 
ciated  with  the  metastable  state  causes  a  very  low  transition 
temperature  /iPX  ~  G(l)  and  gives  rise  to  an  additional 
statistical  factor  exp (SQ)  slowing  down  the  transitions  over 
the  barrier.  This  additional  factor  appeal's  as  a  prefactor  in  the 
Kramers  rate  calculation  [26];  in  the  model  considered  here, 
this  factor  is  exponential  and  needs  to  be  included  to  correctly 
describe  the  scaling  of  the  classical  transition  rate  with  N. 


APPENDIX  B:  EIGENSTATES  OVERLAP 

One  of  the  characteristics  associated  with  complexity  of 
a  given  problem  for  the  quantum  annealing  algorithm  is  the 
overlap  of  the  eigenstate  wave  functions  at  the  beginning 
and  the  end  of  the  algorithm.  We  analyze  it  for  our  model, 
Eqs.  (1)  and  (4).  The  initial  state  s  =  0  is  characterized  by 
the  maximal  x  projection  of  the  total  spin, 


N 

2 


X 


The  overlap  of  this  state  with  the  solution  of  the  classical 
problem,  the  state  fully  polarized  along  z  axis,  is 


X(N/ 2-K\0) 
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